# $Id: rmse.plot.r,v 1.9 2006/04/23 02:44:11 jkatz Exp $



#############################################################
##                                                         ##
##      Plots as T varies                                  ##
##                                                         ##
#############################################################

# Load Data
load("rmse_sim_by_T.Rdata")

### Figure 1: Beta by T
pdf(file="fig1.pdf",width=3.00,height=3.00)
par(las=1,cex=0.5,oma=c(1,1,0,0),pty="m",xaxs="i",yaxs="i")
plot(lowess(rmse.lm[,1],rmse.lm[,2],f=.25),
     axes=FALSE,
     xlim=c(0,50),
     ylim=c(0,2),
     xlab="T",
     ylab="RMSE",
     type="l",
     lty=0)
axis(1,at=seq(0,50,10))
axis(2)
lines(lowess(rmse.lm[,1], rmse.lm[,2],f=.25),  lty=1,lwd=1)
lines(lowess(rmse.lmeML[,1], rmse.lmeML[,2],f=.25),  lty=3,lwd=1)
text(25,1.35,"Pooled OLS",cex=0.75)
text(25,0.30,"RCM (ML)",cex=0.75)
arrows(25,1.30,25,1.05,length=0.03,lwd=0.5)
arrows(25,0.35,25,0.60,length=0.03,lwd=0.5)
dev.off()


### Figure 2: Beta_i by T
pdf(file="fig2.pdf",width=3.00,height=3.00)
par(las=1,cex=0.5,oma=c(1,1,0,0),pty="m",lwd=0.5,xaxs="i",yaxs="i")
plot(rmse.lm[,1],rmse.lm[,4],
     axes=FALSE,
     xlim=c(0,50),
     ylim=c(0,8),
     xlab="T",
     ylab="RMSE",
     type="l",
     lty=0)
axis(1,at=seq(0,50,10))
axis(2)
lines(rmse.lm[,1],rmse.lm[,4],lty=1,lwd=1)
lines(rmse.lmList[,1],rmse.lmList[,4], lty=5,lwd=1)
lines(rmse.lmeML[,1], rmse.lmeML[,4],  lty=3,lwd=1)
text(18,6,"Unit-by-Unit OLS",cex=0.75)
arrows(12,6,7,6,length=0.03,lwd=0.5)
text(20,1.15,"Pooled OLS",cex=0.75)
arrows(20,1.30,20,1.90,length=0.03,lwd=0.5)
text(6,0.25,"RCM (ML)",cex=0.75)
arrows(6,0.40,6,1.00,length=0.03,lwd=0.5)
dev.off()

#############################################################
##                                                         ##
##      Plots as sd(beta_i)                                ##
##                                                         ##
#############################################################


# Load Data
load("rmse_sim_by_sd_beta.Rdata")

### Figure 3: Beta_i by sd(beta_i)
pdf(file="fig3.pdf",width=3.00,height=3.00)
par(las=1,cex=0.5,oma=c(1,1,0,0),pty="m",lwd=0.5,xaxs="i",yaxs="i")
plot(rmse.lm[,1],rmse.lm[,4],
     axes=FALSE,
     xlim=c(0,5),
     ylim=c(0,6),
     xlab=expression(sigma[beta[i]]),
     ylab="RMSE",
     type="l",
     lty=0)
axis(1)
axis(2)
lines(rmse.lm[,1],rmse.lm[,4],lty=1,lwd=1)
lines(rmse.lmList[,1],rmse.lmList[,4], lty=5,lwd=1)
lines(rmse.lmeML[,1], rmse.lmeML[,4],  lty=3,lwd=1)
text(3.50,5.00,"Pooled OLS",cex=0.75)
arrows(3.50,4.75,3.50,4.10,length=0.03,lwd=0.5)
text(1.00,3.45,"Unit-by-Unit OLS",cex=0.75)
arrows(1.00,3.20,1.00,2.55,length=0.03,lwd=0.5)
text(2.50,1.50,"RCM (ML)",cex=0.75)
arrows(2.50,1.25,2.50,0.65,length=0.03,lwd=0.5)
dev.off()

#############################################################
##                                                         ##
##      Plots as population is split                       ##
##                                                         ##
#############################################################
# Load Data
load("rmse_sim_split_pop1.Rdata")

### Figure 4: Beta_i in split population design
pdf(file="fig4.pdf",width=3.00,height=3.00)
par(las=1,cex=0.5,oma=c(1,1,0,0),pty="m",lwd=0.5,xaxs="i",yaxs="i")
plot(rmse.lm[,1],rmse.lm[,4],
     axes=FALSE,
     xlim=c(0,5),
     ylim=c(0,3.5),
     xlab=expression(paste("Size of two outlier ",beta[i])),
     ylab="RMSE",
     type="l",
     lty=0)
axis(1)
axis(2)
lines(rmse.lm[,1],rmse.lm[,4],lty=1,lwd=1)
lines(rmse.lmList[,1],rmse.lmList[,4], lty=5,lwd=1)
lines(rmse.lmeML[,1], rmse.lmeML[,4],  lty=3,lwd=1)
text(1,2.80,"Unit-by-Unit OLS",cex=0.75)
arrows(1,2.70,1,2.50,length=0.03,lwd=0.5)
text(2.50,1.35,"Pooled OLS",cex=0.75)
arrows(2.50,1.25,2.50,1.05,length=0.03,lwd=0.5)
text(4,0.95,"RCM (ML)",cex=0.75)
arrows(4,0.85,4,0.65,length=0.03,lwd=0.5)
dev.off()

#############################################################
##                                                         ##
##      Plots as T varies but with Gamma Distribution      ##
##                                                         ##
#############################################################

# Load Data
load("rmse_sim_by_T_gamma.Rdata")

### Figure 5: Beta by T
pdf(file="fig5.pdf",width=3.00,height=3.00)
par(las=1,cex=0.5,oma=c(1,1,0,0),pty="m",xaxs="i",yaxs="i")
plot(lowess(rmse.lm[,1],rmse.lm[,2],f=.25),
     axes=FALSE,
     xlim=c(0,50),
     ylim=c(0,2.5),
     xlab="T",
     ylab="RMSE",
     type="l",
     lty=0)
axis(1,at=seq(0,50,10))
axis(2)
lines(lowess(rmse.lm[,1], rmse.lm[,2],f=.15),  lty=1,lwd=1)
lines(lowess(rmse.lmeML[,1], rmse.lmeML[,2],f=.15),  lty=3,lwd=1)
text(25,1.35,"Pooled OLS",cex=0.75)
text(25,0.30,"RCM (ML)",cex=0.75)
arrows(25,1.30,25,1.05,length=0.03,lwd=0.5)
arrows(25,0.35,25,0.55,length=0.03,lwd=0.5)
dev.off()


### Figure 6: Beta_i by T
pdf(file="fig6.pdf",width=3.00,height=3.00)
par(las=1,cex=0.5,oma=c(1,1,0,0),pty="m",lwd=0.5,xaxs="i",yaxs="i")
plot(rmse.lm[,1],rmse.lm[,4],
     axes=FALSE,
     xlim=c(0,50),
     ylim=c(0,8),
     xlab="T",
     ylab="RMSE",
     type="l",
     lty=0)
axis(1,at=seq(0,50,10))
axis(2)
lines(rmse.lm[,1],rmse.lm[,4],lty=1,lwd=1)
lines(rmse.lmList[,1],rmse.lmList[,4], lty=5,lwd=1)
lines(rmse.lmeML[,1], rmse.lmeML[,4],  lty=3,lwd=1)
text(18,6,"Unit-by-Unit OLS",cex=0.75)
arrows(11.5,6,7.5,6,length=0.03,lwd=0.5)
text(20,1.15,"Pooled OLS",cex=0.75)
arrows(20,1.30,20,1.90,length=0.03,lwd=0.5)
text(6,0.25,"RCM (ML)",cex=0.75)
arrows(6,0.50,6,0.90,length=0.03,lwd=0.5)

## legend(29,7.90,
##        c("Unit-by-Unit OLS","Pooled OLS","RCM (ML)"),
##        lty=c(5,1,3),lwd=c(1,1,1),cex=0.75)
dev.off()

#############################################################
##                                                         ##
##      Plots for Appendix with FGLS                       ##
##                                                         ##
#############################################################

# Load Data
load("rmse_sim_by_T.Rdata")


### Figure A.1: Beta_i RCM ML vs FGLS
pdf(file="figA1.pdf",width=3.00,height=3.00)
par(las=1,cex=0.50,oma=c(1,1,0,0),pty="m",lwd=0.5,xaxs="i",yaxs="i")
plot(rmse.lmeFGLS[,1],rmse.lmeFGLS[,4],
     axes=FALSE,
     xlim=c(0,50),
     ylim=c(0,4),
     xlab="T",
     ylab="RMSE",
     type="l",
     lty=0)
axis(1,at=seq(0,50,10))
axis(2)
lines(rmse.lmeML[,1],  rmse.lmeML[,4], lty=1,lwd=1)
lines(rmse.lmeFGLS[,1], rmse.lmeFGLS[,4],lty=5,lwd=1)
#lines(rmse.lmeREML[,1],rmse.lmeREML[,4],lty=5)
text(20,2.15,"RCM (FGLS)",cex=0.75)
arrows(20,2.05,20,1.80,length=0.03,lwd=0.5)
text(30,0.85,"RCM (ML)",cex=0.75)
arrows(30,0.75,30,0.50,length=0.03,lwd=0.5)
dev.off()

### Figure A.2: sigma_{beta_i} RCM ML vs FGLS
pdf(file="figA2.pdf",width=3.00,height=3.00)
par(las=1,cex=0.5,oma=c(1,1,0,0),pty="m",xaxs="i",yaxs="i")
plot(rmse.lmeFGLS[,1],rmse.lmeFGLS[,3],
     axes=FALSE,
     xlim=c(0,50),
     ylim=c(0,60),
     xlab="T",
     ylab="RMSE",
     type="l",
     lty=0)
axis(1,at=seq(0,50,10))
axis(2)
lines(rmse.lmeML[,1],  rmse.lmeML[,3]+1, lty=1,lwd=1)
lines(rmse.lmeFGLS[,1], rmse.lmeFGLS[,3],lty=5,lwd=1)
#lines(rmse.lmeREML[,1],rmse.lmeREML[,4],lty=5)
text(20,15,"RCM (FGLS)",cex=0.75)
arrows(20,14,20,9,length=0.03,lwd=0.5)
text(10,8,"RCM (ML)",cex=0.75)
arrows(10,7,10,2,length=0.03,lwd=0.5)

dev.off()




